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Abstract 

An example of full solution of the inverse scattering problem on the half line (from 0 to oo) 
is presented. For this purpose, a simple analytically solvable model system (Morse potential) is 
used, which is expected to be a reasonable approximation to a real potential. First one calculates 
all spectral characteristics for the hxed model system. This way one gets all the necessary input 
data (otherwise unobtainable) to implement powerful methods of the inverse scattering theory. 
In this paper, the multi-step procedure to solve the Marchenko integral equation is described in 
full details. Excellent performance of the method is demonstrated and its combination with the 
Marchenko differential equation is discussed. In addition to the main results, several important 
analytic properties of the Morse potential are unveiled. For example, a simple analytic algorithm 
to calculate the phase shift is derived. 
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I. INTRODUCTION 


Strict mathematical criteria for the unique solution of the one-dimensional inverse scat¬ 
tering problem have been formulated long ago, thanks to the outstanding contributions by 
many researchers^^— (seo^^ for a historical overview). However, in spite of the perfect con¬ 
sistency and mathematical beauty of the concept, these rigorous criteria have rarely been 
used in practice, because the necessary input data are very difficult to obtain. Indeed, one 
has to know the full energy spectrum of the bound states <0 (n = 0,1, 2,..., — 1) and 

the full energy dependence (from 0 to cx)) of the phase shift 5{E) for the scattering states 
> 0). Even if the mentioned obligatory requirements were fulhlled (which is never the 
case), this would only be sufficient to construct an A^-parameter family of phase-equivalent 
and isospectral potentials. One has to somehow £x N additional parameters, the so-called 
norming constants, in order to ascertain the potential uniquely. On the other hand, the 



FIG. 1. An illustration to the physical background of the proposed approach. 


shape of a real potential is not arbitrary, but often it looks like the curve in Fig. [TJ In this 
hgure, a comparison is made between an ah initio potential for Ne 2 molecule!^ and its Morse 
approximant, calculated according to the formula^ 
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taking D = 42.153 K— and Re = 3.0895 A—. The Morse curve in Fig. [T]was calculated 
with ro = {2mD) = 0.2388 A, a = 0.4879/ro, and Ar = r — Re. It has only one bound 

state: Eq/D = —0.5716, and practically the same value has been determined for the single 
discrete level of the Ne 2 molecule, both theoretically^ and experimentally^. 

Fig. [T]was presented for illustrative purposes, to continue the above discussion. Namely, 
in principle, the unknown norming constants for a real potential can be treated as variational 
parameters. One may try to guess their correct values by htting to Eq. ([T]). Moreover, one 
might even think about using this simple model potential to examine different solution 
schemes of the inverse scattering theory, and this is the main idea developed in this paper. 
Morse potential has many unique analytic properties, which makes it especially suitable for 
modelling. In particular, for this model system one can ascertain all the necessary spectral 
characteristics to any desired accuracy. Thereafter, one can use them as input data, instead 
of getting these data from real experiments. 

The idea to invert an inverse problem is, of course, tautological but not at all meaningless. 
Indeed, despite the immense computational power which is now available to the researchers, 
one hardly could hnd any example of complete solution of the inverse scattering problem 
on the half line. The main reason for this dehciency is the lack of the necessary input data, 
and consequently, there seems to be no motivation to perform such studies. Then, what 
is the motivation for the present study? First of all, in author’s opinion there is a debt of 
honor to be paid to the researchers who developed this beautiful theory. Nowadays people 
tend to consider the inverse scattering a ’pedagogical’ problem which is of little scientihc 
interest. Unfortunately, such an ignorant view is not based on the real knowledge, but is a 
mere belief. Implementation of the methods of the inverse scattering theory is not at all a 
trivial task. On the contrary, this is a complex and computationally demanding multi-step 
procedure which has to be performed with utmost accuracy. In this paper, we are going 
to describe this procedure for a very simple model system, to make it understandable for 
a wide audience of potential readers. This can be considered a hrst step towards the real 
implementation of the method. Another, more specihc motive for this work is the obvious 
need to clarify some conceptual details, as will be explained in Sec. II. 

Performing an applicability test of the Marchenko method is not the only issue addressed 
in this paper. In author’s opinion, the model system itself deserves special attention. Useful 
properties of the Morse potential, as well as the solution to the related Schrodinger equation 
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already given by Philip Morse himselfl^, are well known to physicists. Since the potential 
is shape-invariant, this solution can even be found by pure algebraic means, using the Gen- 
denshtein’s recipe^^. However, in this context we are talking about the physically correct 
linear combination of the two special solutions to the Schrodinger equation. There are lots 
of systems which can be analyzed in exactly the same way (the simplest one is the harmonic 
oscillator), but within this pattern the feature which makes Morse potential really unique 
is overlooked. Namely, for a Morse potential the two linearly independent solutions to the 
Schrodinger equation can always be easily ascertained analytically (even ii D < 0!). This is 
not the case for a harmonic oscillator or any other popular model system. 

Consequently, the approach based on Eq. ([T]) is not just an option among many others, 
but it is often the most appropriate choice for modelling a real quantum system. Moreover, 
one may construct a more reliable model consisting of several Morse-type components, still 
preserving the exact solubility of the Schrodinger equation. Such technically rather complex 
developments are not analyzed in this paper (see, e.g.,— and^ for more details), because 
here we concentrate on solving the inverse not the direct Schrodinger problem. For this 
reason, the model has been chosen to be as simple as possible. 

Nevertheless, as a matter of fact, the scattering properties of the Morse potential have 
not found much attention so far, and there still is some work to do. In Sec. Ill, a simple 
analytic formula for calculating the phase shift is derived. Unexpectedly, in addition to this 
specific result, there is a pure mathematical outcome which is directly related to Eq. ©■ 
Namely, a very accurate algorithm for evaluating the Riemann-Siegel function (see^^) was 
obtained, as a kind of bonus for calculating the phase shift. 

The remaining parts of the paper are organized as follows. The details of calculating the 
kernel of the Marchenko equation are described in Sec. IV and the solution method itself is 
under examination in Sec. V. Finally, Sec. VI concludes the work. 


II. MARCHENKO INTEGRAL EQUATION 


The basis for this paper is the integral equation 


CXJ 

A(r, t) = AQ{r + ^) + J ’S)Ho(s -F t)ds, t>r, 


( 2 ) 
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derived by Marchenko^. The Marchenko method is preferred, because the kernel of Eq. 
(1^ is directly related to the main spectral characteristics of the system, while an alterna¬ 
tive approach, the Gelfand-Levitan method^, requires an additional calculation of the dost 
function. 

The kernel of Eq. ([2]) reads 


Ao{x) 


TT / [S{k) - l]exp{ikx)dk 

_ - 

k = \Je/C, 7n = \J-En/C, S{k) = exp [2i6{k )], 


( 3 ) 


where C = h?/ (2m) and Sn are the norming constants for the dost solution of the related 
Schrodinger equation: 


= exp(- 7 „r) -h J A{r,x)exp{--fnx)dx, 

r 

so that —>■ exp(— 7 „r) as r —)■ oo, and 

OO 

4 / [/(^7n,r)]^dr = 1. 

0 

If one is able to solve Eq. ([2]) then the potential can be easily determined: 

,dA{r, r) 
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The key formulas were given in full detail, in order to clarify the issue mentioned above. 
Namely, in several highly appreciated overviews either the Marchenko equation itself or its 


20, p. 73; ll2|, p. 79; l2l|, p. 732; 
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kernel or both of them are given incorrectly (see Refs, 
p. 178). The dispute concerns the sign of the second term in Eq. ((31), which is ’plus’ in 
all these overviews. According to Marchenko’s original works (see, e.g.. Refs. 
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, p. 62; 
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p. 218), on the contrary, there should be ’minus’ as in Eq. ((3]), and this is dehnitely true. 
Indeed, the original Marchenko equation can be rewritten as 

OO 

A{r, t) = -Bo{r + 't) - J ^(a s)Bo{s + t)ds, 
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^o(a;) = — / [I - S{k)]exp{ikx)dk+ ^slexp{-'yrix). 
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Consequently, taking An (x) = —Bq{x), one comes to Eqs. 


sign given in Refs. Il2|, [20 


22 


). Most likely, the wrong 


is just a human error. However, it is misleading (especially in 
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the pedagogical context) and should be corrected. On the other hand, it is quite surprising 
that this error has not even been noticed until now. Indirectly, it means that no one has 
attempted to put the Marchenko method into practice, which provides another motivation 
for the present study. 

Let us specify the problem and the scheme for its solution. First, one calculates the spec¬ 
tral characteristics for a simple model potential described by Eq. ([1]). Then one determines 
the kernel of the Marchenko equation according to Eq. ([3]). Finally, one solves Eq. ([2]) and 
calculates the related potential according to Eq. ([H]). If the procedure is successful then 
the initial potential is recovered. Throughout this paper, dimensionless units for the energy 
and radial coordinate will be used, taking tq = \Jh^/{2mD) = 1, D = 1, and consequently, 
(7 = 1 . 

For a Morse potential, the norming constants can be easily ascertained analytically. 
Let us hx, for example, a = 2/3. Then the system has only one bound state: 

En = -D{1 - = -D{1 - (7) 

a 2a 

where a = D/C/a. Consequently, in our unit system, a = 1/a and 70 = 2/3 = a. 
Therefore, this particular case is especially suitable for illustrative purposes and only this 
case will be examined further on. (The same approach, with minor adjustments concerning 
the different values for a and 70 , can be applied to the model potential shown in Fig. [1]). 
To conclusively hx the model, let us assume that the minimum of the potential is located 
at Re = 2.5ro. 

At hrst sight it may seem that the system has another level exactly at the zero point. 
Indeed, if one takes a = 2/3 then a = 3/2 and, according to Eq. ([7]), Ei = 0. Actually, 
however, this is not the case, because the coordinate here ranges from 0 to 00 , not from 
—00 to 00 as for a classical Morse oscillator. Therefore, as can be easily proved, the possible 
additional level is shifted slightly above (not below) the dissociation limit, becoming a scat¬ 
tering state. Another nuance is that, strictly speaking, one should form a linear combination 
of both solutions of the Schrodinger equation, to fix the correct regular solution. Again, it 
can be shown (see, e.g.,~) that the second special solution is of next to no importance, so 
that the physical solution reads 


To = Noexp{-y/2)y 


( 8 ) 
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as in the ’normal’ Morse case. Here Nq is the related norming constant and a new dimen¬ 
sionless coordinate y = 2a exp [—a (r — Re)] was introduced. 

The parameter Nq as well as the norming constant for the dost solution (0]) can be easily 
determined: 


So = 


l/o« 


yo = 2aexp {aRe) 


(9) 


1 - / exp{-y)ydy 
yo 

The second term in the denominator is negligible, so that, to a high accuracy, Sq = y^a = 
6exp(10/3). This is the value which we are trying to recover. However, to clarify the sign 
problem described above, the quantity Sq will be treated as a variational (not dehnitely 
positive!) parameter. We will see that its expected value is indeed the correct value. 


III. PHASE SHIFT FOR THE MORSE POTENTIAL 


Let us see how to determine the full energy dependence of the phase shift, subjected to 
a strong constraint 

d(0) — d(cx)) = Ntt (10) 

{N = 1 in our case), according to the Levinson theoren>i'^. The analytic background of the 
problem has been described elsewhere^ (see Eqs. (23) to (29) there). There are several 
formally equivalent approaches to building the regular solution of the Schrodinger equation 


for a scattering state. For example, we may proceed from Eq. (25) of Ref. 
some analytical work, we get the following special solutions: 

Ti(/c,r) = 15 {a,i/3;y)\ cos{q:o + arg [S {a,i/3;y)] - kr} , 
T 2 (/c,r) = \S {a, i/3;y)\sin {ao +arg [S {a,i/3;y)] - kr} , 
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FIG. 2. Phase shift for the potential Eq. ([T|) depending on the wave number (parameters are 
specified in Sec. II). Lower graphs demonstrate the asymptotic behavior of 5{k). 


{N = 1 in our case), F being the gamma function. 

From the regularity requirement 'b(O) = y4i'Fi(0) + 742'F2(0) = 0, one gets the ratio of 
the coefficients 


^2 


tan [tto + arg (S'o)] , Sq = S (a, if3;yo). 
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FIG. 3. Wave functions for the scattering states at extremely small energies. The same set of wave 
functions is shown in all graphs and their insets, but note that the coordinate ranges are essentially 
different. 


On the other hand, since S{a, if3] ?/) ^ 1 as r —)■ oo^^, one gets a formula for the phase shift: 




^ costto + sinao 

tan -. 

^ sin do — cos ao ‘ 


(16) 


and consequently, a general expression for the wave function th = Ai^i+A 2^2 (for simplicity 
we drop the arguments of the S'-functions) 

^2 151 


\Ii(/c, r) = 


cos [ao + arg (S'o)] 


sin [kr + arg (S'o) — arg (S)] 
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The last formula is valid for any r, and also in the limit r —)■ oo. From this one immediately 
concludes that A 2 = — cos [ao + arg (So)] and Ai = sin [ao + arg (S'o)] , which means that 

r) = \S\ sin {kr + arg (S'o) — arg (S')} , (17) 

and most importantly, 

5(fc) = arg(S'o). (18) 

This is the simple formula mentioned in Sec. I. Combined with Eq. flT^ . it enables 
to easily ascertain the phase shift for any scattering state and thus, its full dependence on 
energy (or wave number). The result for the specihed Morse potential is shown in Fig. [2] 
(main graph). Note that the /s-scale there is logarithmic, so that the figure involves 18 orders 
of magnitude in the energy space. Full agreement with the Levinson theorem fITU]) can be 
seen as well. 

Fig. [3] demonstrates the behavior of the scattering wave functions as E —)• 0. In general, 
the wave functions can be calculated using Eqs. dnD, dni) and flTSl) . Note, however, that 
the S'-functions are dehned as^ 


S{a, c; x) = exp{—x/2)^{—a + c + 1/2, 2c + 1; x). 


(19) 


where 


, , ax a(a + l)a;^ 

$(a, c;x) = l + — + + .... 


( 20 ) 


l!c 2!c(c + 1) 

Therefore, if E —)■ 0 and consequently, i(3 —)■ 0, then it is more appropriate to explicitly use 
Eq. fl20l) . As a result, one gets the following rapidly converging expression: 


S{a,if3;y) = exp{-y/2) |l -y + if3y 
which was actually used to compose Fig. [3l 


0 ! 


1 ! 


iV 


iW (3!) 
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iV 


( 21 ) 


As the energy is extremely low, linear coordinate dependence can be seen in a wide region, 
which is in full agreement with the well-known rule 

S{k) = Nn — arctan(/cao), (22) 

where oq is the scattering length. From Eqs. flTTll . (|T8ll and (l22il it follows that 

\h(fc, r) —)■ sin [kr + S{k)] —)■ sin (kr) -|- tan [5(A))] —)■ fc (r — oq) , 


and this is indeed the case, while oq = —4312.06224. 

Practically the same scattering length is obtained from the real /c-dependence of the phase 
shift (see the lower left graph in Fig. [2]), using Eqs. flT^ . ffT8|) and ffT^ - fl2^ . 
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A. An algorithm for the Riemann-Siegel function 


Usually, to calculate 6{k) there is no need to make any use of Eqs. In some 


special cases, however, the convergence of the series 
appropriate to combine Eqs. (25) and (28) from Ref. 
analytical work, one comes to the following expression: 


may be slow. Then it may be 
18 l . As a result, again after some 


tan [tto + arg (S')] = 

'( 

where = |(*/d + 1) (*/d + 2)... (i/3 + A)|, 
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(N+iy+y^ 

^ l\y 


[N+lf+y^] [(V+2)2+/32 


(24) 


2!y2 
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and we assumed that A = a — 1/2. 

To exploit Eqs. fl2^ - (l2T|) . one takes y = |/o, uses Eq. (B, and hnally gets the phase 
shift according to Eq. (ITSll . The only problem is how to ascertain the parameter B dehned 
by Eq. flTTl) . Fortunately, we can apply to the theory of gamma functions (see^. Sec. 1.2): 


T{l/2 + z)T{l/2- z) = 


TT 


COS {nz) ’ 


r{ 2 z) = ^^r{z)r(i /2 + z) = 

^ ^ ^ ^ ' ’ cos(7r^) T{l/2-z) 

Taking here 2z = l/2 + i/3, one easily gets the desired result: 


B = arg 


r(z/3) 


= arg [T(i/3 + 1/2)] + 2/3 In 2 = /3 In (Stt) + arctan 


tanh 


'/37r' 


(25) 


+ 29(/3), 
(26) 


where 

0(/3) = arg [r(l/4 + t/3/2)] - /3/2 In vr (27) 

is the Riemann-Siegel theta function, whose behavior is well known as /3 —)■ 0 or /3 —)■ oo. 
Namely^, 


(9 


/33 


/3 ^ 0 : e{(3) = -^ [27 + TT + 2 In (Svr)] + ^ K + 28C(3)J - 
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Here 7 = 0.5772156649... is the Euler-Mascheroni constant and C(3) = 1.2020569032..., (^(5) 
1.0369277551... are the corresponding values of the Riemann zeta function. In practice, Eq. 
fl2S]) can be used if /3 < 0.1, while Eq. fl29]) is recommended for (3 > 20. For the intermediate 
region [3 G [0.1,20] an innovation can be introduced. Namely, combining Eq. with a 
useful formula for calculating arg [r(i/5 -|- 1/2)]^, we get 


20 (/ 3 ) = /3 


In 


(1 + 4 / 32 ) 


47r 


- 1 


— arctan 


tanh 
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wherein 


00 d ^ 

I{I3)= J ^cotht —e“*sin(2/9t)-^ = J e~^ sm{7r^)f{t)dt, 
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coth t 


T = — 
2 /?’ 


— e 


_.^coth(t + T) — _ 2 ^coth(t -I- 2T) — ^ 
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t t + T t + 2T 

The integral /(/3) can be conveniently evaluated numerically or even analytically^: 
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71—1 


f _ 1 \ 71—10271 JD 

HP) = E 4,4 = i:(-i)‘ 


71=1 


(2n)(2u - 1)(1 -l- 4 / 32 ) 2 n- 


k=0 



2k-\-l 


(30) 


(31) 

(32) 


(33) 


Here Bn denotes the n-th order Bernoulli number {Bi = 1/6, B 2 = 1/30, etc.). 

Thus, returning to the main subject of this section, we have found another option for 
calculating the phase shift, based on Eqs. (ITSll . (IT^ . (l23ll . (l26jl and (IHOil - 


6{k) = arctan 


Foaj^l3 exp{yo)y/ 


-{2N+1) 


— tan(ao) 


, Fo = F {N,if3;yo ). 


(34) 


1 — tan (cto) 

This little excursion once again demonstrates wonderful analytic properties of the Morse 
potential, providing an instructive example how purely physical considerations may lead to 
a pure mathematical result. 


IV. SCATTERING PART OF THE KERNEL 

As the necessary spectral data are now available, we can start solution of the inverse 
problem. First step is to perform the Fourier transform to £x the scattering part of the 
kernel ([ 3 ]): 

00 

As{x) = — [ [5'(fc) — 1] exp{ikx)dk = — ^ ( 35 ) 

271 J 71 
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( 36 ) 


f(x) = j sm[26{k)]sm(kx)dk, g{x) =‘^ J sin^ [6{k)]cos{kx)dk. 

0 0 

It can be shown that As{—x) = [f{x) — g{x)] /vr, so that the inverse Fonrier transform 
wonld result in the following formulas: 

OO OD 

1 /* 2 /* 

sin^ [5(fc)] = — g{x) cos{kx)dx, sin[25(fc)] = — f{x) sm{kx)dx. (37) 

71 J 71 J 

0 0 

Calculations according to Eqs. represent the most challenging part of the overall 

procedure, because the functions f{x) and g{x) must be ascertained in a wide range with 
an extreme accuracy. It means, in particular, that the fast Fourier transform techniques are 
useless for our purposes. Fortunately, one can use an asymptotic formulai^^ 

(5(fc) = y +-p +-^ + ..., fc^oo, (38) 




FIG. 4. Demonstration of the excellent accuracy of the Fourier transforms performed according to 
Eqs. (I35])-(I37D. 
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FIG. 5. Coordinate dependence of the function f{x) calculated according to Eq. 


where the coefficients 01 , 03 , 05 ,... are directly related to the potential and its derivatives 

00 

(for example, Oi = —(2G)“^ / V{r)dr). Thus the asymptotic part of the integrals fl5^ can 

0 

be calculated analytically. In addition, one can make use of the Filon’s quadrature formulas 


(see Ref. 


29l. p. 890), which are perfectly suitable if x is large. The Filon’s formulas have 


been used for the range k G [20,100] when x > 20, while for a; < 20 a more common (but 
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FIG. 6. Coordinate dependence of the function g{x) calculated according to Eq. ([36]) • 


very accurate) approach was used. Namely, the domain was divided into sufficiently small 
intervals where 64-point Gauss-Legendre quadrature formula is appropriate to the purpose. 

The described procedure ensures the correct asymptotic behavior of both functions de- 
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FIG. 7. Coordinate dependence of the scattering part of the kernel according to Eq. ([3]) (two 
upper graphs) and the full kernel for the examined Marchenko equation (bottom graph). 

fined in Eq. flSBD . Namely, 

J{x> ~ ^l/oexp(— y) + bexp{-cx), g{x) ^ -yoexp {—^) -6exp(-ca;), (39) 

where b = 7.252681534782e-4, c = 2.315574387346e-4 are some characteristic constants. 
Quite remarkably, these formulas hold to better than 10“^° accuracy for x > 30. The second 
terms in Eq. fl3^ are not important for the following analysis, since they mutually compen¬ 
sate each other. On the other hand, without these additional terms one cannot correctly 
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perform the inverse Fourier transform according to Eq. (1371) . In fact, there is no need for this 
inverse transform as well, but it is recommended, to be sure that the calculated functions 
/(x) and g{x) are reliable. Fig. 0] demonstrates that this is indeed the case. 

FigureslH inland [7] give more detailed information about the solution procedure. (Although 
the curves in these hgures may look like solid lines, they are assemblies of dots which 
represent the actual results of calculations). Looking at Figs. IMZl one can understand why 
the utmost accuracy is needed in performing the Fourier transform. This is most clearly 
seen in the lowest graph of Fig. [7] (the green dashed line there shows A(0)). Indeed, if 
the Fourier transform is not sufficiently correct then the important information about the 
scattering properties of the system (shown in hgures [5l [6] and in the two upper graphs of 
Fig. [7]) may easily be lost against the background of the bound states’ contribution to the 
kernel, which dominates at small values of the argument. Thus the whole procedure would 
become meaningless. 

As the result of the Fourier transform according to Eq. (I36|) . one gets two smooth 
functions that can be very accurately (with absolute error less than 10“®) approximated by 
piecewise rational functions 

{ /(x) a + 6x + cx^ + dx^ + ex^ + /x® 

1 + gx + hx"^ + ix^ + jx* + fcx® + /x® ’ 

whose parameters are given in Tables [I] and [TTl Using these parameters and applying Eq. 
for larger arguments, one can reproduce all curves shown in Figs. IMZl 

V. SOLUTION OF MARCHENKO EQUATION 

Having completed the preliminary work, there remains the last step - the actual solution 
of the Marchenko equation ([2]). To some surprise, this is technically much easier than 
determining the kernel (0]). For any hxed r, let us define T(x) = A{r, r + x), t = r + x, and 
s = r + y. Eq. ([2]) can then be rewritten as 

OO 

T(x) = Ao(2r + x) + y Ao(2r + x + y)T{y)dy k, Ao(2r + x) + ^ WiF{Xi), (40) 

0 * 

where Fiij) = Aq (2r + x + y)T {y), while Xj and Wi are the nodes and weights of an 
appropriate quadrature formula, respectively. 
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TABLE I. The fitting parameters for the curves shown in Fig. [SI 


Parameters of the rational fit for the function f{x) 



0 < X < 0.2 

0.2 < X < 2.5 

2.5 < X < 10 

10 < X < 30 

a 

-8.0742834752365 

-8.07451927770804 0.710177479533 

-1.318246783662 

b 

-139.70663281 

21.09396832161 

0.078971528944 

0.350978767195 

c 

-1964.01791531 

-20.103272203664 

-0.276506389451 

-0.025997992558 

d 

9295.18176575 

10.475159850219 

0.065047583073 

8.813724334413e-4 

e 

-8951.34532192 

-2.595277468302 

-4.695098041353e-3 

-1.492795622311e-5 

f 

0 

0.196451244097 

1.138463611195e-4 

1.117794694102e-7 

9 

19.59774821 

-0.316168857071 

-0.711060782617 

-0.345296699762 

h 

287.27716444 

0.804977555407 

0.277492349329 

0.068087218353 

i 

-511.76899191 

-0.412291982417 

-0.053313022152 

-6.848766446371e-3 

j 

-357.91228265 

0.355627589898 

6.620633459645e-3 

4.849331544015e-4 

k 

-747.01469965 

-0.117132459825 

-3.875251671984e-4 

-1.655209945751e-5 

1 

0 

0.024943017758 

1.561945583680e-5 

4.455108455684e-7 


Eq. fHOj) has been solved using the Nystrom method (see, e.g., Ref. yO, P- 782), which 
is perfectly suitable for our purposes. As the asymptotic behavior of the kernel is known, a 
good idea is to split the integral into two parts and discretize them differently. Namely, for 

R 

h = f F{y)dy the 64-point Gauss-Legendre quadrature formula can be successfully used. 
0 

More specihcally, the family of isospectral potentials in Fig. [H] has been calculated with 
R = 15, and splitting the range [0, R] into two subdomains, so that there are 128 nodes 
within this region. 

OO 

To discretize the asymptotic integral R = J F{y)dy, a useful trick is to change the 

R 

variable^: y = R-|-Q!o(1—'w)/(1+m), where ao > 0 is an arbitrary constant. Correspondingly, 
the domain of integration reduces to [—1,1], where one can again use the 64-point Gauss- 
Legendre quadrature formula. Thus 


Jo — 2an 


-1 


F{u) 
(m + 1)' 


,du = Y,w[F{x[), 


(41) 


where x[ = R + ao(l — Xi)/{1 + Xj), w[ = 2aQWi/ (x* -t-1)^, and are the nodes and 

weights of the Gauss-Legendre formula, respectively. 
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TABLE II. The fitting parameters for the curves shown in Fig. O 


a 

b 

c 

d 

e 

f 

g 

h 

i 

j 

k 

I 


Parameters of the rational fit for the function g{x) 

0<x< 0.25 0.25 < X < 2.42 2.42 < x < 9.2 9.2 < x < 42 


3.072175507081882 

19665382.109603 

-57958653.686259 

37126643.692319 

14117338.25322 

-11218702.142292 

6400798.980337 

2752650.386055 

3804036.247907 

-49419.924115 

-268374.263668 

1638562.625682 


3.070830178732 

-0.08759364482 

-24.726764443843 

34.416513599217 

-13.367864496127 

1.423672500917 

3.338766375967 

0.570620534536 

1.953451792127 

-0.164502127028 

0.34471652522 

0.011395750684 


-4.755191843506e-5 

0.149350091364 

-0.015988794423 

-0.032070876143 

6.826423264373e-3 

-2.754362667700e-4 

-1.199750111544 

0.712576159696 

-0.216463815475 

0.03998655984 

-3.895572554791e-3 

2.046893098608e-4 


-0.350846307261 

0.087431813456 

-5.295698261860e-3 

1.228492194728e-4 

-7.497417661209e-7 

-9.379725943814e-9 

-0.403835495838 

0.082504174267 

-8.976469204382e-3 

6.288995558866e-4 

-2.318297714188e-5 

5.388478464471e-7 


The computational procedure itself is simple. Namely, using the discretized version of Eq. 
(HOI) at the specihed node points Xn, one gets a system of linear algebraic equations for the 
quantities T„ = T{xn)- Such a system can be easily solved. In this paper, the Householder 
reduction algorithm's was used for this purpose, because it is numerically more stable than 
Gaussian elimination. Having found the solution, one uses Eq. fHOj) once again, taking x = 0 
to get T(0) = H(r, r). Finally, one calculates the potential according to Eq. (j6]). 

The results of solving the Marchenko equation are shown in Figs. [8]and[9l The curves 
there were calculated with oq = A ■ (1 -f Xo)/(l — Xq) and A = 10000 (this parameter 
characterizes the width of the domain). As can be seen, discretization of Marchenko equation 
with only 192 nodes related to the relevant quadrature formulas is sufficient to accurately 
solve the inverse problem. In addition, one can convince himself/herself that the sign of the 
second term in Eq. ([3]) is indeed ’minus’, while Sq = yla = 6exp(10/3), as expected. 

Fig. [H] demonstrates the dependence of the solution A(r, r + x) on both arguments, r and 
X. In this case, the exact theoretical value for the norming constant (sq = |/o«) was used to 
hx the kernel. Actually, to calculate the potential according to Eq. ([6]), one only needs the 
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FIG. 8. 3D depiction of a solution to Eq. m (see the explanations in the text). 


solution at X = 0, which is shown as a background contour (dash-dot line). 

In Fig. [9l a family of isospectral potentials can be seen. It was obtained by varying the 
norming constant sq in a wide range (from 0 to 100), while the scattering part of the kernel 
remains exactly the same for all these curves. Looking at Fig. [HI one can imagine that 
if the real norming constant is not known, a criterion of ’reasonableness’ could be used to 
§x it. For example, assuming that the potential has only one extremum point (minimum) 
and only one inflection point (as the Morse potential), the correct norming constant can be 
determined with reasonable accuracy. 
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4 6 8 
Dimensionless radial coordinate 


FIG. 9. A family of phase-equivalent isospectral potentials calculated according to Eqs. m- 
(j6]) . The black solid line shows the initial Morse potential, yellow circles correspond to the exact 
theoretical norming constant Sq = y^a = 168.18975, black dashed line indicates the position of the 
energy level Eq = — 7 q = 4/9, and green dotted line demonstrates the result of ’removing’ this level 
(so = 0). 


VI. CONCLUSION 

In this paper, a detailed description of the full solution procedure of the Marchenko 
integral equation was given. To this end, the complete set of the necessary input data 
was used, which was obtained by accurately solving the direct Schrodinger problem for a 
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model system. Such a tautological combination of the direct and inverse scattering problems 
may seem useless for practical purposes, but one has to bear in mind that otherwise the 
full set of necessary data cannot be obtained at all. Unfortunately, the methods of the 
inverse scattering theory are still not implemented as planned. These methods have mainly 
been used to construct families of isospectral potentials^i^, and they have many interesting 
applications in the soliton theory (see, e.g.,— but their capabilities have not been fully 
exploited. Hopefully, the results of this work clearly demonstrate the real power of the 
Marchenko method, albeit the necessary input data were obtained computationally, not 
experimentally. On the other hand, the model system used for this purpose, the well-known 
Morse potential, continues to surprise us by its undiscovered analytical properties. For 
example, in this paper a simple formula for calculating the phase shift was derived, along 
with the new analytic algorithm for calculating the Riemann-Siegel function, which is a pure 
mathematical result. 


How could the results of this work be useful for further studies? An interesting option 
is to combine Eq. (Ej) with the Marchenko differential equation (see, e.g.. Ref. Il2|, p. 78), 
assuming that the real potential can be expressed as a sum V{r) = Vo(^) + AU, where Vo(r) 
is a known model potential whose spectral characteristics are known as well. Thus according 
to Eq. ([6]), Vo(r) = —2 [Ao(r,r)]^, where A^ir^t) is the corresponding solution to Eq. (|2]). 
Then, assuming that V{r) = —2 [A(r, r)]^ and A{r,t) = Ao{r,t) + AA, one gets an equation 


dr"^ 



AA(r, t) 


Vo{r)AA{r, t) + AV (r)Ao(u t). 


(42) 


If V^ir) is sufficiently close to the real potential, the right side of Eq. (H2ll would contain 
two small quantities, AU and AA, which, in principle, can be determined self-consistently 
along with solving this equation. 
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